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Using empirical potential molecular dynamics we compute 
dynamical matrix eigenvalues and eigenvectors for a 4096 
atom model of amorphous silicon and a set of models with 
voids of different size based on it. This information is then em- 
ployed to study the localization properties of the low-energy 
vibrational states, calculate the specific heat C(T) and ex- 
amine the low-temperature properties of our models usually 
attributed to the presence of tunneling states in amorphous 
silicon. The results of our calculations for C(T) and "excess 
specific heat bulge" in the C(T)/T S vs. T graph for void- 
less a-Si appear to be in good agreement with experiment; 
moreover our investigation shows that the presence of local- 
ized low-energy excitations in the vibrational spectrum of our 
models with voids strongly manifests itself as a sharp peak in 
C(T)/T 3 dependence at T < 3K. To our knowledge this is 
the first numerical simulation that provides adequate agree- 
ment with experiment for the very low-temperature proper- 
ties of specific heat in disordered systems within the limits of 
harmonic approximation. 



I. INTRODUCTION 

The main goal of the computational project presented 
below was to construct a set of realistic models for de- 
vice quality a-Si (a-Si:H) material, containing nanovoids 
in the structure, and explore the vibrational properties of 
these models, especially the localization patterns of the 
low-energy states that emerge after introducing a void 
into silicon continuous random network (CRN). Mod- 
ern computational facilities allow us to study models of 
up to several hundreds of atoms with ab initio methods 
and thousands of atoms if we switch to empirical tech- 
niques, which is especially helpful for realistic large scale 
modeling of amorphous materials. We can perform sim- 
ulated quenching and annealing for these models, look- 
ing for their most energetically favorable geometry, then 
we can calculate the dynamical matrices for our systems 
and diagonalize them exactly, receiving their eigenval- 
ues together with the conjugate eigenvectors. This data 
gives us the ability to construct the vibrational density 
of states (VDOS) for a given model as well as to look at 
the localization properties of individual eigenstates of its 
dynamical matrix. The calculation of specific heat can be 
seen as a natural extension of these techniques because 
we can do it with a little effort once we have obtained 



the VDOS information for the model. 

In Section II of this paper we present the details of our 
model construction scheme for a large (based on 4096 
atom model) family of models for a-Si with voids and 
introduce a set of methods we employ for geometry op- 
timization, dynamical matrix and specific heat calcula- 
tions. In Section III we discuss the results of our cal- 
culations of vibrational properties and specific heat for 
the models considered and finally in Section IV we sum- 
marize our findings about the influence of localized low- 
energy vibrational modes on the thermodynamical prop- 
erties of amorphous silicon. 



II. MODEL CONSTRUCTION AND 
COMPUTATIONAL PROCEDURES 

We use the Djordjevic, Thorpe and Wooten 4096 atom 
model 1 for a-Si (referred to as DTW in what follows) 
constructed with the Wooten, Winer and Weaire bond 
switching algorithm 2 as a base for building a family of 
models with voids. The length of the side of a cubic 
supercell for this model is approximately 43 A. As our 
first step we optimize the geometry of the basic model 
which results only in minor network rearrangements; this 
relaxed geometry is then used for producing all of the 
following models with voids. To cut out a void we pick 
an arbitrary atom in the network and remove it as well 
as the consecutive spherical shells of its neighbors. We 
find that our results do not depend much on which atom 
we select for this procedure. 

By applying this procedure we have built three models 
with voids of different diameter: a 4091 atom model with 
a "small void" (only one atom and four of its nearest 
neighbors removed) — a void of approximately 5 A in 
diameter, 4069 atom "medium void" model with 10 A 
void and 4008 atom "large void" model with 15 A void. 
We refrained from building models with even larger voids 
to prevent possible interaction of a void with its own 
"ghost" images in the neighboring periodically translated 
supercells. 

Every model with voids was then quenched to minimize 
the forces acting on atoms in the network. The atomic 
forces must be small for the application of the harmonic 
approximation for the total energy of the system, which 
is required for the dynamical matrix calculation. After 
the dynamical matrix calculation for every model men- 
tioned above, the eigenvalue and eigenvector data was 
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used to produce VDOS and inverse participation ratio 3 
(IPR) graphs for the model, calculate its specific heat and 
visualize the spatial localization/delocalization charac- 
teristics of some of its vibrational modes (the ones which 
behavior we found most interesting). 

For geometry optimization of the models (simulated 
quenching) and dynamical matrix calculation we em- 
ploy a molecular dynamics code "Estrelle" developed by 
the authors of this paper. The code is built around an 
empirical environment-dependent interatomic potential 
(EDIP), which has been recently introduced by Bazant 
and Kaxiras 4-6 . In general, this potential inherits the 
well established Stillinger- Weber 7 format for two- and 
three-body interactions, but now these interaction terms 
depend on the local atomic environment through an ef- 
fective coordination parameter. Our testing results for 
EDIP and its performance in comparison to our previous 
ab initio calculations 8 are described elsewhere 9 . 

The force tolerance threshold in simulated quenching 
mode of our MD program is set to be 0.01 cV/A for 
all its applications we discuss here. The dynamical ma- 
trix for any given model is computed by displacing every 
atom in the cell by 0.03 A in three orthogonal directions 
and calculating the originating forces on all the atoms 
in the system 10 . For a system of thousands atoms the 
dynamical matrix is very large which, under normal cir- 
cumstances, causes problems in storing it on disk or in 
computer memory. Fortunately the dynamical matrix 
is also very sparse, because in most cases the displace- 
ment of a single atom generates significant forces only 
on its closest neighbors but not in the whole supercell. 
We extensively exploit this localization of dynamical ma- 
trices in our calculations, discarding terms smaller than 
10~ 4 eVA" a.u.m. , which is a good compromise be- 
tween accuracy and compactness of the output. Once 
the sparse dynamical matrix for the system is obtained 
we use a separate routine to exactly diagonalize the whole 
matrix and obtain all of the eigenvalues and eigenvectors. 
Again, for the same reasons as already mentioned above, 
we do not write out all of the eigenvectors (however, we 
do keep all their IPRs) but rather only those that ex- 
hibit properties we look for: (i) small energy/eigenvalue 
(less than 200 cm" 1 ) and (ii) relatively high IPR, show- 
ing that the vibrational mode we are dealing with is lo- 
calized. 

To create VDOS graphs we employ a Gaussian rep- 
resentation for S(E — Ei), where Ei,i = 1 . . . N are the 
eigenvalues of the dynamical matrix and N = 3N atom s- 
The width of broadening is 20 cm" 1 for the full scale 
graphs and 0.1 cm" 1 for the close-ups of the low-energy 
region. The vibrational activity colormaps for the "low- 
energy, high IPR" modes are prepared in same way that 
has been already described in Sec. II of Ref. 8: the set of 
individual atomic IPRs is computed and then every atom 
is assigned a certain color according to its displacement 
from the equilibrium position. 

It is relatively easy to obtain C{T) dependence for the 



model if the VDOS information for it is available 11 : 
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where VDOS g(E) is normalized to unity. Nevertheless 
one thing should be treated with caution: the model 
VDOS one usually has is relevant for a system of finite 
size (i. c. our supercell). Vibrational excitations with 
wavelengths longer than the size of the supercell cannot 
be excited in this model and are consequently missing in 
its VDOS data. In order to receive the precise values for 
C(T) one should correct VDOS for the infinite size of the 
system. In our case it is done in the following fashion: 
all the delocalized (acoustic) vibrational modes of energy 
less than 20 cm -1 are cut out and substituted by a weak 
parabolic tail aE 2 in the routine to compute the VDOS. 
Parameter a can be obtained from a calculation of the 
elastic constants of the model 12 but in this investigation 
we use a more simple approach, fitting a to provide a 
smooth transition between the low-energy parabolic tail 
and the rest of VDOS. 



III. DISCUSSION OF RESULTS 
A. Vibrational properties and localization 

We begin this section by presenting the results of our 
calculations of the low-energy regions 13 for VDOS and 
IPR for all the four models introduced above, shown in 
Fig. 1. For the sake of simplicity we do not present the 
colormaps for vibrational modes for 4096 atom family of 
models in this paper, however a set of colormaps for the 
most interesting localized excitations in these models is 
available for download over the World Wide Web 14 . 

From Fig. 1 we can see that the large models for a- 
Si (both, with and without voids) exhibit quite a com- 
plicated vibrational behavior, much more complex than 
that of smaller 216 atom based families of models, we 
have studied before 8 ' 9 . The most important difference 
here is that a-Si model without voids has two localized 
low-energy modes that are associated with strained re- 
gions of silicon network: we have checked bond lengths 
and bond angles for the atoms in these regions and found 
that the modes generally localize on atoms with bond an- 
gles deviating from the perfect tetrahedral angle by more 
than 30 degrees. 

Consequently, now we have two types of phonon traps 
in our models with voids — the voids themselves and 
strained regions of the network. Keeping this in mind we 
can attempt to introduce a rough classification of the lo- 
calized low-energy modes according to the type of phonon 
trap they fall into. First, we can see a significant number 
of vibrational modes in our models with voids that gen- 
erally show the same kind of localization properties that 
we have reported earlier 8 : they are exponentially local- 
ized with the center of localization positioned to the side 
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of the void. We classify these excitations as void type 
modes. Secondly, the modes we might attribute to the 
strained network region phonon trap type in models with 
voids exhibit a different kind of behavior in comparison to 
the voidless model. These modes do not localize exactly 
on the strained regions in the supercell; instead they form 
a string extended between one of these strained regions 
and the void 15 . The possible explanation of this behavior 
is that these modes can be regarded as a superposition 
of void type modes and the localized excitations in the 
model without voids. In our classification we call them 
mixed type modes. We have to stress once again that 
the classification we propose is only approximate and is 
based mostly on the colormaps (i. e. pictures) we get for 
our models not on rigorous mathematical arguments. We 
must also add that all the low-energy modes, that appear 
localized in our finite models, will be pseudolocalized in 
an infinite sample 16 . 



with voids modes of different types dominate in the low- 
energy region. In the "small void" model a mode with the 
highest IPR at 10.58 cm -1 is of void type, but the suc- 
ceeding three modes with high IPR at 14.43, 18.25 and 
20.97 cm -1 are of strongly pronounced mixed type. In 
the "medium void" model, to the contrary, all three low- 
energy localized modes at 5.89, 6.12 and 8.13 cm -1 are 
of mixed type. The mode with strong void type behavior 
is also present but it is shifted to 17.97 cm -1 . Finally in 
the "large void" model modes at 2.34 and 6.10 cm' 1 arc 
of void type and all the others, including a strongly local- 
ized mode at 10.28 cm' 1 , exhibit mixed type behavior. 
We speculate that the network strain and geometrical pe- 
culiarities of any given model play a more important role 
in shaping the energy and type distribution of its local- 
ized vibrational modes than the actual size of the void 
— at least for the models with voids of comparable sizes, 
as we have here. 
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FIG. 1. Low-energy VDOS and IPR snapshots for 4096 
atom DTW model without voids (upper left set of panels), 
4091 atom "small void" model (upper right), 4069 atom 
"medium void" model (lower left) and 4008 atom "large void" 
model (lower right). 

We must admit that in our current investigation we 
were not able to find any simple connection between the 
size of the void and the energy and type of resulting lo- 
calized modes. Our data shows that for different models 



B. Specific heat 

In this section we present our results for the calcula- 
tions of specific heat C{T) for the family of 4096 atom 
models. The overall temperature dependence for specific 
heat for all of our models is in good agreement with Du- 
long and Petit's law at high temperatures and Debye's 
law at low temperatures; our calculation also produces 
approximately the correct Debye temperature for a-Si. 
For the room temperature (300K) we receive practically 
the same value for specific heat for all of our models: 
19.7 JK^mor 1 . 

In the left panel of Fig. 2 the C(T)/T 3 low-temperature 
dependence for our models is presented. The most strik- 
ing feature in this graph is the presence of sharp peaks 
at T < 3K in the curves for the models containing 
voids. The model without voids does not have this peak, 
although it does demonstrate the presence of the well 
known excess specific heat bulge, the position and height 
of which are in qualitative agreement with experiment 17 
as well as with recent computational results of Feldman, 
Allen and Bickham 16 . All of our models with voids also 
have the excess specific heat bulges at approximately the 
same position, but comparing to the low-tcmperature 
peaks their intensities are about an order of magnitude 
smaller. We were not able to find any experimental data 
for specific heat measurements in a- Si at temperatures 
below 2K, but in order to make some general comparison 
to experiment for these new low-temperature features we 
obtain (which should be generic for any disordered sys- 
tem containing voids), we provide the experimental curve 
for vitreous silica 18 in our graph. 

Unlike the previous void size vs. energy situation, we 
can find a clear connection between the presence of low- 
energy localized modes in vibrational spectrum of the 
model and the height (or even absence) of the low- 
temperature peak in its C(T)/T 3 graph. For the "small 
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void" model we have a localized mode at 1.19 cmT x (here 
and below, see Fig. 1) — the lowest energy at which we 
can see localized excitations in all our models — and 
the highest peak in C(T)/T 3 dependence. The "large 
void" model has its lowest energy localized excitation at 
2.34 cm -1 and the peak of smaller height comparing to 
the previous model. The "medium void" model has two 
localized states but only at approximately 6 cm -1 and 
peak that is even less pronounced than in case of the 
first two models. Finally the model without voids has no 
localized states with energy lower than 11 cm -1 and no 
low-temperature peak whatsoever. 
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FIG. 2. Left panel: the low-temperature C(T)/T 3 depen- 
dence for the DTW models with and without voids. The 
experimental curve for vitreous silica is taken from Ref. 18. 
Right panel: the curves for 4091 atom "small void" model 
before and after clipping of the lowest energy localized mode 
eigenvalue. Curve for the model without voids is also shown 
for reference. 

In order to investigate this connection in more detail 
we have performed a simple numerical experiment, which 
results are shown in the right panel of Fig. 2. We have 
clipped the eigenvalue at 1.19cm -1 from the eigenvalue 
set for the "small void" model and recalculated its VDOS 
and C(T) receiving no low-temperature peak in C(T)/T 3 
graph, much like in the situation with the model with- 
out voids. In our opinion these results provide enough 
evidence to attribute the existence of low-temperature 
(T < 3K) peak in C(T)/T 3 dependence for the model to 
the presence of localized low-energy (E ~ 1 — 6 cm ) 
vibrational excitations — in our case produced by voids 
in its spectrum. 

Finally we must note that the localized vibrational ex- 
citations we see, although having rather low energies, 
are not tunneling states 18 , that are nonharmonic by na- 
ture and can not be obtained in harmonic approximation 
calculation. We do not claim that the whole tunneling 
states theory is incorrect, we rather propose an alter- 
native mechanism that explains the same experimental 
data. It seems that any mechanism that creates ad- 
ditional density of vibrational states (be this tunneling 
states or low-energy localized "void" vibrations in porous 



materials) at very low energies will produce the same ef- 
fect on low-temperature specific heat behavior. In or- 
der to find out which mechanism of the two mentioned 
above actually works in real material, an experimental 
investigation of low-temperature thermal properties and 
simultaneously geometrical quality (i.e. presence of de- 
fects, voids, strained regions) of this material should be 
carried out. The works of X. Liu et al. 19,20 or Coeck and 
Laermans 21 for amorphous silicon can be regarded as the 
closest examples here. 

IV. CONCLUSIONS 

We have studied vibrational and thermodynamical 
properties of 4096 atom DTW model for amorphous 
silicon and the family of models with voids based on 
it, employing Bazant-Kaxiras environment-dependent in- 
teratomic potential and empirical MD technique. We 
have found that the models with voids posses a complex 
spectrum of localized low-energy excitations that can be 
roughly divided into two groups — void and mixed type 
modes — according to their localization patterns. Our 
calculations show that there is no simple connection be- 
tween the size of the void and the energy and type of 
its localized modes. It is most probable that not only 
the size of the void but also its local geometrical envi- 
ronment as well as strain distribution in the neighboring 
regions of the network play a paramount role in shaping 
the low-energy vibrational spectrum of the system. We 
have constructed specific heat C{T) plots for our models, 
that appear to be in good agreement with experiment. 
We have also plotted out our models' C(T)/T 3 depen- 
dences for the low-temperature region, which seem to be 
in adequate agreement with experimental and other com- 
putational results for T > 3K (the excess specific heat 
bulge) and predict new interesting features, undoubtedly 
connected with vibrational properties of voids present in 
the system, at lower temperatures. We must stress that 
our results are correct for model materials with a uniform 
distribution of voids of one and the same size, which is 
of course impossible to produce in real material. Never- 
theless, employing our model data we can predict that in 
real material the localized low-energy vibrational states, 
connected to voids of different sizes, will fill out a band 
which will alter the parabolic VDOS tail properties at 
small energies and consequently manifest itself by chang- 
ing the specific heat C(T)/T 3 dependence. 
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